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^ ; Abstract 

,, An algorithm for the evaluation of the complex exponential function is proposed 

which is quasi-linear in time and linear in space. This algorithm is based on a 
Q 1 modified binary splitting method for the hypergeometric series and a modified Karat- 

c/2 ■ suba method for the fast evaluation of the exponential function. The time complexity 

' of this algorithm is equal to that of the ordinary algorithm for the evaluation of the 

exponential function based on the series expansion: 0(M(n) log(n) 2 ). 

>: 

CN . 1. Introduction. In this paper we introduce a measure of the space complexity of 

calculations on a Schonhage machine [T] and give an upper bound for the time and space 
fSJ 1 complexity of a proposed algorithm for the computation of the exponential function of 

qq ■ a complex argument in each area \z\ < 2 P where p is a natural number, p > 0, on the 

<0 Schonhage machine. 

The Schonhage machine is in fact an ordinary computer. Therefore, we describe an 
• • ■ algorithm which is quasi-linear in time and linear in space on an ordinary computer. 

. . Hence, the phrase 'fast algorithm' refers to evaluations on a Schonhage machine. In 

particular, the estimate 0(n log (n) log log (n)) refers to this machine [I]. 

Basic information on dyadic rational numbers and constructive real numbers and 
functions can be found in [2j. The notation Sch(FQLINTIME / /LINSPACE) will 
be used for the class of algorithms which are quasi-linear in time and linear in space on 
a Schonhage machine. Quasi-linear means that the complexity function is bounded by 
0(nlog(n) fc ) for some k. 

From now on, n will denote the length of the record of accuracy 2~ n of dyadic 
rational approximations; x will be used for a real argument, z wiil be used for a complex 
argument. We will use log (A;) for logarithms base 2. 

The basic subject of interest is the algorithms for the evaluation of elementary 
functions based on series expansions, as such algorithms are important for practical 
computer science due to the relative simplicity of their implementation. 

Algorithms for the fast evaluation of the exponential function and some other ele- 
mentary functions with a time complexity of 0(M{n) log(n) 2 ) are offered in [3| (where 
M(n) denotes the complexity of multiplication of n-bit integers); the space used by 
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these algorithms is bounded by 0(nlog(n)), as will be shown below. In [4j, algorithms 
for the calculation of elementary functions, based on Taylor series, which are quasi- 
linear in time are considered; the space used by the algorithms from [3] is bounded by 
a quasi-linear function as the classical binary splitting method uses 0(nlog(n)) space 
for the intermediate results. 

That is, elementary functions are computable using quasi-linear time and quasi- 
linear space. There is a question: whether it is possible to evaluate elementary functions 
using quasi-linear in time and linear in space algorithms based on series expansions? 

This paper shows that the answer to this question for the complex exponential 
function and some other complex elementary functions is Yes. We use the combi- 
nation of two algorithms for the construction of an algorithm of complexity class 
Sch(FQLINTIME / /LINSPACE) for the evaluation of the exponential function: 
a modified binary splitting method for the evaluation of the hypergeometric series, and 
a modified Karatsuba method [3] for the fast evaluation of the exponential function. 

Note that the residual sum of the series ([7]) satisfies the following inequality: 

\R u (r)\ < C2-( rlos M +m ); 

here r = m2~ u+1 . Therefore, |i?jy(r)| = 0(2 -mlog ( m )) doesn't hold, and we cannot get 
the result about FLINSPACE computability of exp(x) from this equation. 

2. Description of the computation model (machine Schonhage). This 
machine, introduced in [1], operates on symbols of the alphabet S = {0, 1, . . . , 2 A_1 } 
and sequences of such symbols (we can take, for example, a constant A equal to 32). 
The machine consists of arrays Tq, . . . ,T T to read and write symbols from S, registers 
A, B, C, M for arithmetic operations, and a control unit CPU. The arrays are infinite 
in both directions. For each array there is a pointer pi to the current symbol written 
in the array. The record < p + j > means the symbol referenced by pointer p + j. 
There is also an additional register Y which is a pointer to the current operation, and 
an optional array S which serves as the stack of recursive calls; S is infinite in one 
direction. A bit register E acts as an overflow register for arithmetic operations. 

A program for Schonhage consists of several modules written in the language TPAL, 
which is similar to an assembly language for a RISC processor. In TPAL there are 
commands for loading a symbol written in an array into a register, for reading a symbol 
from a register and writing it to an array, for increasing and decreasing the content of 
a register, the shift command, the call and return from a procedure commands, the 
jump to a label command, and the conditional jump command. Integers on which the 
machine operates are encoded as symbol sequences in the alphabet E: 

a = a + ai2 A + a 2 (2 A ) 2 + . . . + a fc _i(2 A ) fc - 1 , < a. L < 2 A - 1. 

The sign bit is written in the symbol which is before the senior symbol afc-i- 

Schonhage can call procedures and perform recursive calls. After the return from a 
procedure the memory occupied by the parameters and local variables is released. 

The time computational complexity of an algorithm on Schonhage is defined as the 
number of instructions in the language TPAL. Arithmetic operations on symbols and 
calls of procedures are counted as a constant number of steps. The memory used in an 
array during the calculation is defined as the maximum of the number of array elements 
involved in the calculation. 

As constructive functions are functions that compute approximations of functions 
using approximations of arguments, we define the oracle machine Schonhage. This 
machine has some oracle functions that compute approximations of arguments; the 
machine calculates approximations of a function using these approximations of argu- 
ments. A request to an oracle is written in array Tq as the record of an accuracy of 
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the computation; approximations of arguments are recorded in array To too. A query 
to an oracle is treated as one operation in the time computational complexity of the 
oracle machine Schonhage. 

Definition 1. The space computational complexity of an algorithm on the oracle ma- 
chine Schonhage is defined as the sum of the memory used for all the arrays plus the 
maximum of the memory used for the stack. 

3. Constructive complex numbers and functions. A complex number z' 
such that \z — z'\ < 2~ n is called an approximation of the complex number z with 
accuracy 2~ n . 



Suppose that there are complex numbers oj = x + iy, u' = x' + iy' such that 
x - x'\ < 2~( n+1 ) and \y - y'\ < 2~( n+1 ). Then 



That is, to calculate an approximation of complex number uj with accuracy 2 n it is 
sufficient to calculate an approximation of the real and imaginary parts of the complex 



We say that a sequence (f> : N — > D X D, where (j){n) = (<p x (n),4> y (n)), D is the 
set of dyadic rational numbers, converges dyadic-rationally to the complex number z 
if for any n 6 N the following holds: prec(4> x (n)) = n + 2, prec(4> y (n)) = n + 2, and 
\z(n)—z\ < 2~ n , where z(n) = (f) x (n)+i4> y (n). The set of all functions 4> which converge 
dyadic-rationally to a number z is denoted by CF Z . A complex number z is called a 
CF constructive complex number if CF Z contains a computable function (p. 

Definition 2. A complex number z € C is called a Sch(FQLINTIME / /LINSPACE) 

constructive complex number if there exists a function (f) € CF Z which belongs to the 
class Sch(FQLINTIME / /LINSPACE) . 

Let / be a function f(z) : A — > C, where A = {z € C : \z\ < R} is an area in the 
set of complex numbers. 

Definition 3. The function f(z) is called a Sch(FQLINTIME//LINSPACE) con- 
structive complex function in the area A if for any z from this area there is a function 
tp from CF f(z) which belongs to the class Sch(FQLINTIME / /LINSPACE) . 

Note that to calculate the values of a constructive complex function of the argument 
z = x + iy we need to specify functions u(x,y) = Re{f{z)) and v(x,y) = Im(f(z)), 
and in order to calculate these functions we need to have two oracle functions that 
correspond to the real and imaginary parts of the argument. 

4. Binary splitting method. This method is used to calculate the values of 
series with rational coefficients, in particular, to calculate the hypergeometric series of 
the form 



where a, b, p, and q are polynomials with integer coefficients. Linearly convergent 
hypergeometric series are used to calculate many constants of analysis and elementary 
functions at rational points; this series is linearly convergent if its partial sum 



oj — a/ 1 = 



v 7 ^ 



x') 2 + {y- y') 2 < V2 ■ 2-2(™+ 1 ) < 2- n . 



(l) 



number with accuracy 2 ('■ 





(2) 
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where fJ.(k) is a linear function of k, differs from the exact value by not more than 2 



S - S(n(k))\ < 2 



In its classical variant, the binary splitting method works as follows. Put k\ = fJ>(k). 
We consider the partial sum ([2]) for some integers i\ and i 2 , < i% < ki, < i 2 < fci, 
k < i 2 - 



We calculate P(h,i 2 ) = p(h) ■ • -K^), Q(h,h) = q(h) ■ ■ ■ B{i x ,h) = b(i x ) . . . b(i 2 ), 
and T(ix,i 2 ) = B(ix,i 2 )Q(ii,i 2 )S(ii,i 2 ). If i\ = i 2 then these values are calculated di- 
rectly. Otherwise, the series is divided into two parts, left and right, and P(i x ,i 2 ), 
Qihik), an d B(i\,i 2 ) are calculated for each part recursively. Then the values ob- 
tained are combined: 



The algorithm starts with i\ = 0, i 2 = k\. After calculating T(0, ki), B(0,k\), and 
Q(0, ki), we divide T(0, k\) by B(0, k\)Q(0, k\) to get the result with the given accuracy. 
We write out the binary splitting method explicitly 

Algorithm BmSplit. 

Approximate value of the partial sum ([2]) with accuracy 2 fe . 

Input: Record of accuracy 2~ k . 

Output: Approximate value of l^fy with accuracy 2~ fe . 
Description: 

1) ki := fj,(k); 

2) [P,Q,B,T] := BinSplitRecurs(0,ki); 

3) perform division r := -j^ with accuracy 2~ fc ; 

4) return result r. 

This algorithm uses the following subalgorithm computing recursively the values P, 
Q, B, and T. 

Algorithm BinSplitRecurs. 
Calculation of P, Q, B, and T. 

Input: Bounds i\, i 2 of the interval. 

Output: Tuple [P(h,i 2 ), Q(h,i 2 ), B(i x ,i 2 ), T(i 1 ,i 2 )]. 

Description: 

1) if ii = i 2 then T := a{ii)p{i\) and return tuple \p(ii), q{h), b(ii), T]; 



3) calculate [P[,Qi,Bi,Ti] := B in Split Recur s^ijimid); 

4) calculate [P r , Q r , B r , T r ] := B in Split Recur s(i m id,i 2 )', 

5) T := B r Q r Ti + BiP{T r and return tuple [P t P r , QiQ r , B t B r , T). 

The lengths of T(0, k\) and B(0, ki)Q(0, k\) are proportional to fclog(/c); there- 
fore the binary splitting method is quasi-linear in space; the time complexity of this 
algorithm is 0(M(k) log(k) 2 ) @]. 




P(h,i 2 ) = PiP r , Q(h,i 2 ) = QiQr, B(ix,i 2 ) = BiB. 
T(i 1 ,i 2 ) = B r Q r T l + B l P l T r . 



(3) 



_ h+i2. 
2 ' 
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5. Karatsuba's method for fast evaluation of exp(x). We consider the 
Taylor series of the real exponential function 



00 i 1 

X „ X X X" 



^ = E-r = 1 + TT + ^ + --- + ^ + --- ( 4 ) 

at xq, —\ + 2~ m < Xo < \ — 2~ m . We compute the value of this series with accuracy 
2 _n3 , 1x3 > 7, at the dyadic rational point x m , \x m — xq\ < 2~ m , — -? < x m < -?. Let k 
be the smallest value such that 

n 3 + 1 < 2 k , m = 2 k+1 . (5) 

We represent x m = ±0.000304 . . . a m a m +i as 

X m = ±0.000304 + ±0.000005060708 + • • • + ±0.00 . . . 0a m _ 2 fc +1 a m „ 2 ft +2 • • • a mO-m+l = 

ft , ft ^ ft , , ft+i 

= ^4+^8+^16+--- + ~2m~ = 72 ± 73 ± • • • + 7fc+l> 

where ft = ±a 3 o 4 , ft = ±a 5 a 6 o 7 a 8 , . .., ft+i = ±a m _ 2 k +1 a m _ 2 k +2 ■ ■ ■ a m a m+ i; 
7i/ = ft2 -2 , 2 < 1/ < A;+l; ft/ is a 2^~ 1 -digit number. We write exp(x m ) as a product: 

exp(a; m ) = exp(7 2 ) exp(7 3 ) . . . exp(7 fc+ i). (6) 

In Karatsuba's method (FEE method, fast evaluation of the exponent; this method is 
also known as the Brent trick [5]) the values exp(7„) are then calculated using a Taylor 
series 

exp( 7 ,) = 1 + yj^ + ^7 + • • • + ^7 + Ru{r) = £„ + R„(r); (7) 
here r = m2~ u+1 . As for the residual sum, the following inequality is satisfied [3J: 

|ft| r+1 

|i? ^ r)l < 2 (r + l)r2(''+ 1 ) 2 ^ 
then |i?„(r)| < 2" m . The values £ ^ are obtained from the formulas 

f - — a -£ b b - r'2 r2 " 

K 

where the integers a v are computed using a sequential process of grouping members of 
the series (JT]). Note that the formula for b v contains the factorial of r, and r varies from 
ttt-2" 1 to m2~ k . Hence there are values which are proportional to 713!, and therefore the 
length of the intermediate results is proportional to n3log(ri3). 



6. Modification of the binary splitting method. We will modify the bi- 
nary splitting method for the evaluation of the hypergeometric series ([7]) so that the 
calculations are in the class Sch(FQLINTIME//LINSPACE). 

We take r = m2~ u+2 1 k\ = log(r), and n = [^-] {k\ is a natural number), and 
write the partial sum of ([7]) as 

P{lv) =<?!+ T 2 [(J2 + T 3 [<7 3 + . . . + T fcl _i[cT fel _i + T fel <7 fc J]], (8) 
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where 

a l = 1 + TT77o77 + „,„o.o„ + ■■■ 



T2 



V.2 2v 2!2 2 ' 2 " (n - l)!2( r i- 1 )- 2 " 



ri !2n-2^' 



a 2 = 1 + 7 , -, No9 „ + 7 ; — 777 , n\^.^ + • • • + 



(n + l)2 2 " (r 1 + l)(r 1 + 2)2 2 - 2 ' / - " ( n + 1) . . . (2 n - l^i- 1 )' 2 " 
ft 1 



(ri + 1) . . . 2r x 2 r i 



CT 3 = l + 77; — ; 777^ + 77; — ; 7777: — , nXn9 . 9 ^ + • • • + 



(2ri + l)2 2 " (2n + l)(2n +2)2 2 - 2 " (2ri + 1) . . . (3ri - l^fa- 1 )- 2 " 



The at are calculated by the classical binary splitting method for the sum ([2]), where 

ji(k) = n - 1, a(i) = 1, fe(i) = 1, 

'l i = (9) 




We estimate the computational complexity of the calculation of at and Tj in the 
following two lemmas. 

Lemma 1. The time complexity of the binary splitting method for the calculation of 
at on the Schonhage machine is bounded above by (rlog(r) + m) log(m) log log(m); the 
space complexity is bounded above by 0(m), where m is given by @. 

Proof. We consider an arbitrary maximal chain of recursive calls derived from the 
calculations in accordance with algorithm BinSplitRecurs. We consider pairs (P,i), 
(Q,i), (B,i), and (T, i) where i is the number of an element of the chain of recursive 
calls and P, Q, B, and T are in the ith element of the chain. We suppose that the 
numbering in the chain begins with its deepest element: i = 1, ... ,q, where q is the 
length of the chain, q < [log(2ri)]. 

Using mathematical induction on i, we show that the length of the representation 
of T in the pair (T,i) satisfies l(T) < 2%u) + 2\ where u = r2 2 ". We note that 
l(Pv) < l(u) since (3 U is a 2^ 1 -digit integer. Next, l[P) < 2%u), 1{Q) < 2%u) for pairs 
(P, i), (Q,i), since increasing i leads to doubling the length; B is always 1. 

The induction begins with i = 1: l(T) < l(u) < 2 1 l(u) + 2 1 . This inequality follows 
from ^ and the formula T = a(i\)p(i\) for (T, 1). Next, the induction step is 

1{T) < 2%u) + 2%u) + 2 l + 1 < 2 i+1 l{u) + 2 i+1 . 

This inequality follows from ([3]): there are two multiplications of numbers of length 
2 l l(u) and 2 l l(u) + 2* and one addition. 

We note that, based on this inequality, T of the pair (T, q) has length 0(m) because 



l(T) < d2H(u) < C 2 -—_(log(r) + 2 U ) < C 2 
log(r) 



m2~ u+l 2 v 

r -\ 

log(r) 



< C 3 m. 



We estimate the time complexity of the calculation of at ■ We must take into account the 
property of the complexity of integer multiplication: 2M(2~ 1 m) < M(m) (quasi-linear 
and polynomial functions satisfy this property). Since at a tree node of recursive calls 
at level i there are C4 multiplications of 2 ?_l numbers of length at most 2 l l{u) + 2 l and 



6 



h < Clog(m), we obtain the following estimate for the number of operations required 
to calculate at- 

Time(a t ) < d ^ 2^ l M(2 i l(u) + 2 l ) < C 5 2<- l M(2 i+1 l(u)) 

i=l i=l 

i=i 

< C 7 qM(2H(u)) < C 8 log{r)M (r + 



i=l i=l 

m 
log(r) 

(here we use the inequality for 2 ? l(u) from the estimate of l(T) for the pair (T, The 
final division gives 0(M(m)) operations. If we use the Schonhage-Strassen algorithm 
for integer multiplication, then 

(in \ 
r + - — — log(r)log(m)loglog(m) = 
log(r)/ (10) 

= C§{r log(r) + m) log(m) log log(m). 

Now we estimate the space complexity of the computation of at- At an element 
of a chain of recursive calls with number i, the amount of memory consumed for the 
temporary variables is Cio(2 l l(u) + 2*). Hence, we conclude that the amount of memory 
in all the simultaneously existing recursive calls in the chain is estimated as follows: 

Space(a t ) < ^ do (2*1 (u) + 2*) < C u (2H(u) + 2 C ) < C 12 m = 0(m). 
i=l 

□ 

Lemma 2. The time complexity of the calculation of Tt using binary splitting method 
on a Schonhage machine is bounded above by (r log(r) + m) log(m) loglog(m); the space 
complexity is bounded above by 0(m) where m is given by ([5]). 

Proof. The estimates of the computational complexity of Tj are the same as those for 
at due to the fact that the inequalities in the proof of Lemma [T] are also suitable for 
Tt (we can calculate the numerator and denominator of at using the binary splitting 
method for products). □ 

We calculate approximate values P{^ u )* with accuracy 2~( m+1 ) by ([8]) using the 
following iterative process: 

/ii(mi) = a* kl , 

hi(mi) = Ofo-i+i + 7fc 1 _ i+2 ^_i, i = l,...,h, (11) 
hi{m\) = hi(mi) + Si\ 

for i = k% we set P{^ p )* = /i^(mi). Here m\ > m (mi will be chosen later), and a* and 
t* are approximations of dj and Tj with accuracy 2~ mi . The hi{mi) are obtained by 
discarding bits q mi+ iq mi+ 2 - - - q mi +j 01 numbers hi(mi) after the binary point starting 
with the (mi + l)th bit: 

\ei\ = \hi(m\) -hi(mi)\ = 0.0. . . 0q mi+ iq mi+2 ■ ■ ■ q mi +j, (12) 

and the sign of £$ is the same as the sign of hi(m\) (it is clear that |ej| < 2~ mi ). 

We note that the accuracy of exp(7„) is 2" m since we compute £* = P(^ u )* with 
accuracy 2~( m+1 ) and 

| exp(7,)* - exp( 7 ,)| < \£ - £„| + \R„(r)\ < 2~( m+1 ) + 2~( m+1 ) = 2~ m . 
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Lemma 3. For every i G 1 . . . k± 

\h i (m 1 )\<2. (13) 
Proof. We apply mathematical induction on j for hj(m\) using the estimates 

4 1 

| ai| < exp(7„) < -, Tj < 7/ 1 < -. 

The induction base for j is j = 1: |/ti(mi)| < cr^ + 2~ mi < 2. The induction step for 
(j + 1) > 2: 



2 + 2~ mi < 2. 



4 

£ 3 + 



4 +Z 



□ 



Lemma 4. TTie error of the calculation o/fet, (mi) using the scheme (jlip is estimated 
to be 

Proof. We put 

Hi = a kl , Hj = cr fcl + r fel _j +2 i/j_i, 77(1, mi) = |/i»(mi) - 
We use mathematical induction for rj(j,mx) on j. The induction base for j is 1: 

r/(l, mi) = |^i (mi) - Hi I = |< -a kl \< 2~ mi+1 . 
The induction step is (j + 1) > 2: 

r/(j + l,mi) = + r^^+afyCmi) + £ i+ i - o- kl _ {j+1)+1 - T klHj+1)+2 Hj\ 

< \T*hj(mi) - T v hj(mi) + T v hj{mi) - r v Hj\ + 2 ■ 2~ mi 

< 2-" 11 /i j (mi) + 2- 2 r/(j,mi) + 2 • 2~ mi . 

Since (|13p . \hj(mi)\ < 2. By the induction hypothesis, mi) < 2~ mi+J , and so we 
get 

n{j + l,mi) < 2 • 2~ mi + 2- 2 2- mi+ -? + 2 • 2-™ 1 < 2- mi +^' +1 ). 
From A(A;i,mi) = n(ki,mi) we now obtain the required inequality. □ 

Lemma U] implies that it is sufficient to take mi = 2m + 1 to compute P(ju) with 
an accuracy of 2 _ ( m+1 ). 

We denote the algorithm for the calculation of the hypergeometric series using 
scheme ([TT]) by RLinSpaceBinSplit (linear space binary splitting). 

Algorithm RLinSpaceBinSplit. 

The approximate value of the hypergeometric series. 

Input: Record of the accuracy 2~ m . 

Output: The approximate value of (J7J) with accuracy 2~ m . 
Description: 

1) mi := 2m + 1; 
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2) h := 0£ (using the classical binary splitting method with accuracy 2~ mi ); 

3) make a loop through % from 2 to k\\ 

a) calculate V\ := cr^ with accuracy 2~ mi using the classical binary splitting 
method and v 2 := T k 1 -i+2 w ith accuracy 2 _mi , 

b) calculate h := V\ + v 2 h, 

c) assign value h to h rounded in accordance with (|12p ; 

4) write /i on exit. 

We estimate the time computational complexity of this algorithm on a Schonhage 
machine, taking into account that m dependends linearly on 773: 

• log(r) computations of at give the following (from inequality (jlOp ): 

log(m) 

Time(all(at)) < Og(r log(r) + m) log(m) loglog(m) < 

< Ci 3 (mlog(m) 2 loglog(m)) = 0(M(n 3 ) log(n 3 )); 

• 0(log(ri3)) computations of give 0{M(n^) log(ns)); 

• 0(log(n3)) multiplications of numbers of the length 0(n 3 ) give 0(M(n 3 ) log(n 3 )); 

in total we obtain 0{M{n^)\og{n^)). The space complexity of the modified binary 
splitiing method RLinSpaceBinSplit is 0(713) since in all calculations in this algorithm 
we process numbers of length 0(713). 

Proposition 1. The modified binary splitting algorithm RLinSpaceBinSplit belongs to 
the class Sch(FQLINTIME / /LINSPACE) . 

7. Modification of FEE. We construct a modification of the method FEE so 
that our calculations are in class Sch(FQLINTIME // 'LIN SPACE). 

In the classical algorithm FEE the values exp(x m )* are calculated by §6§ using 
a pairwise summation of the numbers exp(7j)*; in FEE we need to keep in memory 
O (log (log (713))) numbers of length 0(713) and, as already mentioned, n 3 ! values are 
processed in this algorithm, so it does not have linear space complexity. 

We calculate exp(x m )* with accuracy 2~ n ' A using the following iterative process: 

h 2 (m) = exp(7 2 )*, 

hi(m) = /ij_i(m)exp(7j)*, i = 2,...,k + l, (14) 
hi(m) = hi(m) + ef, 

for i = k + 1 we put exp(x m )* = hf.j r i{m). Here exp(7j)* are approximations for the 
values exp(7j) with accuracy 2~ m obtained from formula ([7]). The values hi(m) are 
obtained by discarding bits g m +ig m +2 • • • q m +t of the numbers hi(m) after the binary 
point starting with the (m + l)th bit, i.e., 

\si\ = \hi(m) - hi(m)\ =0.0... 0q m+ iq m+2 ■ ■ ■ q m +t, (15) 

and the sign of £j is the same as the sign of hi(m) (it is clear that |£j| < 2~ m ). 
Lemma 5. For every i £ 2 . . . k + 1 

\hi(m)\ < 2 1 - 1 . (16) 
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Proof. We apply mathematical induction on j for hj(m). The induction base is j = 2: 
\h 2 (m)\ < | e X p( 72 )| + 2- m < | + 2- 16 < 2 2 - 1 

(here we take into account that | 7 j| <\, m> 16). The induction step is (j + 1) > 3: 

\h j+1 (m)\ = |/ ij (m)exp( 7i+ i)* + e j+1 \ < %(m)(exp( 7j+ i) + 2~ m ) + 2~ m 
"3 



< 2^ 



-i 



_ + 2 -m 

2 



+ 2" m < 2 J . 



□ 



Lemma 6. The error of the calculation of hk+\{m) using scheme (| 14|) is estimated to 
be 

A(Jfe + l,m) < 2- m+2 ( fc+1 ). 

Proof. We put 

F 2 = exp(7 2 ), Hi = exp( 72 ) exp( 73 ) . . . exp( 7 j), r)(i,m) = \hi(m) - Hi\. 
We use mathematical induction for rj(j,m) on j. The induction base is j = 2: 
77(2, m) = \h 2 (m)-H 2 \ = | exp( 72 )* - exp( 72 )| < 2- m < 2 - m+2 ( 1+1 ). 
The induction step is (j + 1) > 3: 

V(j + l,m) = |/ij(m)exp( 7i+ i)* + e,- + i - -fr, exp( 7i+ i)| 

< |/tj(m)exp( 7 j + i)* - /ij(m)exp( 7j - + i)+ 

ft-j(m) exp( 7j+ i) - Hj exp( 7j+ i)| + 2~ m 

< hj(m)2~ m + 7?(i, m) exp( 7j+1 ) + 2~ m . 

Since (I16p . [/ij(m)| < 2 3 . Then by the induction hypothesis, 7](j,m) < 2~ m+2j , and 
so we get 

r/(j + l,m) < 2 J '- 1 2- m + [2- m+2j '] ~ < 2- m+2 ^ +1 \ 

From A(k + 1, m) = r/(k + 1, m), we now get the required inequality. □ 

Lemma [6] implies that the accuracy 2~ m of the calculation of exp( 7 j)* is sufficient 
to calculate exp(x m )* with accuracy 2 -713 , since 

- m + 2(k + 1) < -n 3 =>• 2 k+l - 2(k + 1) > n 3 and 
2 k+1 - 2(k + 1) > 2 fc+1 2" 1 = 2 k > n 3 + 1 

(here we recall that n 3 + 1 < 2 k ). 

We denote the algorithm of the calculation of the real exponential function using 
scheme (fH| by RLinSpaceFEE (linear space fast exponential evaluation). 
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Algorithm RLinSpaceFEE. 

The approximate value of the real exponential function on the interval 

L 8 ' 8 J 

Input: Record of the accuracy 2~ ns . 

Output: The approximate value exp(x) with accuracy 2~ nz . 

Oracles: (j) x . 

Description: 

1) h := exp(72)* (using algorithm RLinSpaceBinSplit with accuracy 2~ m ); 

2) make a loop through i from 3 to k + 1: 

a) calculate v\ := exp(ji)* using algorithm RLinSpaceBinSplit with accuracy 
2~ m , 

b) calculate h := h ■ v\, 

c) assign the value h to h rounded in accordance with (fT5l) ; 

3) write h on exit. 

The time complexity of this algorithm on a Schonhage is 0{M{n^)\og{n^) 2 ) as 
the algorithm of the calculation of the hypergeometric series RLinSpaceBinSplit uses 
0(M(n^) log(ri3)) operations, and in scheme (|14|) there are 0(log(?73)) such calculations 
and O (log (77.3)) multiplications of numbers of the length 0(77.3); the space complexity 
of RLinSpaceFEE is 0(723) since in all the calculations in this algorithm numbers of 
length 0(77,3) are used. 

Proposition 2. The modified FEE algorithm RLinSpaceFEE for the calculation of the 
exponential function belongs to the class Sch(FQLINTIME / /LINSPACE) . 

9. Calculation of the real function exp(x). Let p be a positive integer, p > 0. 
We compute the function exp(x) with accuracy 2~ n in the interval [— 2 P ,2 P ]. 

We perform the multiplicative reduction of the interval of the complex argument. 
Namely, we take an integer s = 2 P+3 and x' = -; then \x'\ = i.e., x' is in the interval 
[— 2~ 3 < x 1 < 2 -3 ]. Thus the calculation of exp(x) is reduced to the computation of 
exp(x') and then we raise this value to the power s to get the result. It is easy to see 
that the dependency function for m of the accuracy of exp(x') is m = L(n) + C(p), 
where L(n) is a linear function of 77, and C(p) is a constant which is independent of p 
(constant in the sense that it doesn't depend on n). 

Put ?77 > 77i + 3. We then have: x m = x + 9\2~ m , \9\\ < 1; |x | < 2 -3 , and 

|exp(x ) -exp(x m )| = | exp(x ) - exp(x ) exp(#i2~ m )| = exp(x )|l - exp(6 1 2~ m )\. 

Since exp(6»i2- m ) < jz^r, |exp(6»i2- m ) - 1| < < 2~ m+1 . Therefore we have 

the estimate 

|exp(x ) - exp(x m )| < 2 • 2-' m+1 = 2- m+2 < 2~^ +1 \ 

which shows provided the accuracy of the calculation of x m is better than 2~( ni+3 ), 
then one achieves an accuracy 2 - ™ 2 , 772 = rii + 1 for exp(xo). Since m = 2 fc+1 , this 
condition is satisfied. 

Now we need to keep in mind that we calculate the approximate value of exp(rc m )*. 
If this approximation is calculated with accuracy 2 _n3 , 773 = n\ + 1, then 

I exp(x ) - exp(x m )*| < | exp(x ) - exp(x m )| + | exp(x m ) - exp(x m )*| 
^ 2~( n i+ 1 ) _)_ 2 _ ( ni + 1 ) = 2~ ni 

This implies that we can take 7x3 = 77i + 1 in algorithm RLinSpaceFEE. 
We are now ready to describe the algorithm. 
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Algorithm RLinSpaceExp Value. 

The approximate value of the complex exponential function. 

Input: Record of the accuracy 2 _n . 

Output: The approximate value exp(z) with accuracy 2 _n . 
Parameters: <j) x for the argument x. 
Oracles: Constant p. 
Description: 

1) m :=L(n) + C(p); 

2) n 3 := m + 1; 

3) calculate fe, m so that ([5]) holds; 

4) pi := p + 3 

5) s := 2Pi; 

6) compute x* := x (max(l,m — pi)); 

7) perform the reduction of the interval: (x*)' = ^- (the accuracy of the arguments 
will be 2~ w ); 

8) using algorithm RLinSpaceFEE, calculate v := exp(x*)* with accuracy 2~™ 3 ; 

9) write complex number v s to the output. 

The properties of algorithms TLinSpaceBinSplit and TLinSpaceFEE allow us assert the 
following propositions. 

Proposition 3. Algorithm RLinSpaceExpValue of the calculation of the complex expo- 
nential function belongs to the class Sch(FQLINTIME / /LINSPACE) . 

The estimates of the computational complexity of algorithm RLinSpaceExp Value on 
the Schonhage machine are the same as those for algorithm RLinSpaceFEE: that is, 
the time complexity is 0(M{n^) log(ri3) 2 ) and the space complexity is 0(123). If we use 
the Schonhage-Strassen algorithm for integer multiplication, then the time complexity 
of algorithm RLinSpaceExpValue is bounded above by 0(n^ log (n^) 3 log log (^3)). 

Proposition 4. The real function exp(x) is a Sch(FQLINTIME / /LINSPACE) 

constructive real function in any interval [— 2 P ,2 P ]. 

8. Calculation of the function exp(i ■ y). It is easy to show that all the 
algorithms, lemmas, and estimates can be formulated for the evaluation of the function 
exp(i • y). 

1. Calculate ([2]), where p(j) = i ■ Preal(j)] P, Q, B are integers, and T is complex. 
The estimates of the computational complexity are the same as the estimates in 
Lemmas [1] and [2l The hi(m) are obtained by discarding bits q mi +2Qmi+3 ■ ■ ■ Qmi+j 
of the numbers hi(m\) after the binary point starting with the (m\ + 2)th bit. 
Lemmas [3] and ?? are true for the new scheme. Algorithm CLinSpaceBinSplit is 
the same as RLinSpaceFEE. 

2. In the FEE method, we calculate 

exp(i • x m ) = exp(i • 72) exp(i • 73) . . . exp(i • 7^+1). 

The estimate for |i?„(r)| is the same as that for exp(x), that is, the series Eq : 
RealExp : GammaNuSeries converges linearly for exp(i • y). 

3. In the formulas for cjj and Tj, we use i-/3 v and we use algorithm CLinSpaceBinSplit 
for the computation of P(i ■ 7^)*. 
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4. In the scheme for the computation of P(i ■ 7^)*, the hi{m) are obtained by dis- 
carding bits <j , mi+2<?mi+3 • • ■ Qmi+j °f the numbers hi{m\) after the binary point 
starting with the {m\ + 2)th bit; Lemmas and [S] are true for the new scheme; 
algorithm CLinSpaceFEE is the same as RLinSpaceFEE. 

Denote the algorithm for the computation of exp(i ■ y) by CLinSpaceExpValue. 

9. Calculation of the complex function exp(z). Let p be a positive integer, 
p > 0. We compute the function exp(z) with accuracy 2~ n in the area \z\ < 2 P (we 
have \x\ < 2P, |y| < 2P). 

We perform the multiplicative reduction of the interval of the complex argument. 
Namely, we take an integer s = 2 P+3 and z' = |; then |z'| = i.e., z' is in the area 
\z'\ < 2 -3 and both |x| < 2 -3 and |y| < 2 -3 . Thus the calculation of exp(z) is reduced 
to the computation of exp(z') and then we raise this value to the power s to get the 
result. It is easy to see that the dependency function for n\ of the accuracy of exp(z') 
is m = L[n) + C(p), where L(n) is a linear function of n, and C(p) is a constant which 
is independent of p (constant in the sense that it doesn't depend on n). 

Next we consider the function exp(z) in the area \z\ < 2~ 3 . Put Q x = exp(x) and 
C, y = exp(i • y). According to (pQ), we need to calculate ( x and £ y with accuracies of 
2~ n ' 2 , ri2 = n\ + 1 respectively, in order to calculate exp(z) with an accuracy of 2~ ni + 1 . 

Put m > m + 3. We have the following: x m = xq + #i2 -m , |^x| < 1; at that 
|xq| < 2 -3 and 

|exp(x ) -exp(x m )| = |exp(x ) - exp(x ) exp(6»i2" m )| = exp(x )|l - exp(6»i2~ m )|. 

Since exp((9i2- m ) < j^r, |exp(6li2- m ) - 1\ < < 2- m+1 , and therefore we 

have the estimate 

|exp(x ) -exp(x m )| < 2 • 2-' m+1 = 2~ m+2 < 2~^ +1 \ 

which shows that if accuracy of the calculation of x m is better than 2 - ( ni+3 ), then 
one achieves an accuracy of 2~( ni+1 ) for exp(xo)- Since m = 2 fc+1 , this condition is 
satisfied. A similar estimate can be obtained for | exp(i • yo) — exp(i • y m )\. 

Now we need to keep in mind that we calculate the approximate value of exp(z m )*. 
If this approximation is calculated with accuracy 2 _n3 , 723 = n\ + 1, then 

[ exp(z ) - exp(z m )*| < [ exp(z ) - exp(z m )| + | exp(z m ) - exp(z m )*| 
^ 2 _ ( n i+ 1 ) _i_ 2 _ ( n i+ 1 ) = 2 _ni 

This implies that we can take 77,3 = n\ + 1 in algorithms RLinSpaceFEE and CLinSpace- 
FEE. 

We are now ready to describe the basic algorithm. 
Algorithm LinSpaceExp Value. 

The approximate value of the complex exponential function. 

Input: Record of the accuracy 2~ n . 

Output: The approximate value exp(z) with accuracy 2 _n . 
Parameters: (f) x and <j) y for the argument z = x + iy. 
Oracles: Constant p. 
Description: 

1) ni := L(n) + C(p); 

2) n 3 := ni + 1; 
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3) calculate k, m so that ([5]) holds; 

4) pi := p + 3 

5) s := 2 Pl ; 

6) compute x* := </) x (max(l, m — Pi)), y* := J/ (max(l, m — pi)); 

7) perform the reduction of the interval: (x*)' = — , (y*)' = ^- (the accuracy of the 
arguments will be 2 _m ); 

8) using algorithm RLinSpaceFEE, calculate v% := Q* with accuracy 2 _ ™ 3 ; using 
algorithm CLinSpaceFEE, calculate Vi := £*; here the arguments are (x*)' , (y*)'; 

9) write to the output the complex number {y\ + i ■ V2Y ■ 

The properties of algorithms RLinSpaceBinSplit, RLinSpaceFEE, CLinSpaceBinSplit, 
and CLinSpaceFEE allow us to assert the following propositions. 

Proposition 5. Algorithm CLinSpaceExpValue of the calculation of the complex expo- 
nential function belongs to the class Sch(FQLINTIME / /LINSPACE) . 

The estimates of the computational complexity of algorithm CLinSpaceExpValue 
on a Schonhage machine are the same as those of algorithms RLinSpaceFEE and 
CLinSpaceFEE; that is, the time complexity is 0(M(n 3 ) log(ns) 2 ) and the space com- 
plexity is 0(713). If we use the Schonhage-Strassen algorithm for integer multiplica- 
tion, then the time complexity of algorithm CLinSpaceExpValue is bounded above by 
0(n 3 log(n 3 ) 3 log log(n 3 )). 

Theorem 1. The complex function exp(z) is a Sch(FQLINTIME / /LINSPACE) 

constructive complex function in any area \z\ < 2 P . 

10. Computation of the complex functions sin(z), cos(z), sh(z), ch(z). 
Based on the formulas for the trigonometric functions 

sin(z) = - — ^ — = f _^ — , cos(z) = 6 + ^ — , 
we obtain the following: 

Proposition 6. The complex function sin(z) is a Sch(FQLINTIME / /LINSPACE) 

constructive complex function in any area \z\ < 2 P . 

Proposition 7. The complex function cos(z) is a Sch(FQLINTIME / /LINSPACE) 

constructive complex function in any area \z\ < 2 P . 

The following two propositions also follow directly from the formulas for the hyper- 
bolic sine and cosine: 

g * j g * 

Mz) = o , ch(x) = . 



Proposition 8. The complex function sh(z) is a Sch(FQLINTIME / /LINSPACE) 

constructive complex function in any area \z\ < 2 P . 

Proposition 9. The complex function ch(z) is a Sch(FQLINTIME / /LINSPACE) 

constructive complex function in any area \z\ < 2 P . 
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11. Conclusion. Constructed algorithm CLinSpaceExpValue can be used in com- 
puter science as the basis of the Sch(FQLINTIME//LINSPACE) constructive com- 
plex functions exp(z), sin(z), cos(z), sh(z), ch(z), defined on the set of Sch(FQLINTI- 
ME/ /LINSPACE) constructive complex numbers. 

Note also that if we use a simple recursive method for integer multiplication with 
time complexity 0(n log ^), then the time complexity of algorithm CLinSpaceExpValue 
is 0(n lo s( 3 ) log(n) 2 ). 

As future research plans, we could note the problem of the construction of algo- 
rithms based on series expansions for the Sch(FQLINTIME//LINSPACE) com- 
putable analogues of other elementary functions as well as the importance of the prob- 
ably more difficult problem of the construction of computable analogues of elementary 
functions (also based on algorithms that use series expansions) with a time complexity 
of 0(nlog(n) fc ), k < 3, and with linear space complexity. 
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